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(D ' 

QThis paper is concerned with conditions for the existence of oscillations in gene 
I regulatory networks with negative cyclic feedback, where time delays in transcription, 

. translation and translocation process are explicitly considered. The primary goal of 

CN| ' this paper is to propose systematic analysis tools that are useful for a broad class of 

cyclic gene regulatory networks, and to provide novel biological insights. To this end, 
' ■ we adopt a simplified model that is suitable for capturing the essence of a large class of 

' gene regulatory networks. It is first shown that local instability of the unique equilib- 

rium state results in oscillations based on a Poincare-Bendixson type theorem. Then, 
(/3 ' a graphical existence condition, which is equivalent to the local instability of a unique 

O I equilibrium, is derived. Based on the graphical condition, the existence condition is 

analytically presented in terms of biochemical parameters. This allows us to find the 
dimensionless parameters that primarily affect the existence of oscillations, and to pro- 
^ , vide biological insights. The analytic conditions and biological insights are illustrated 

1/^ ■ with two existing biochemical networks, Repressilator and the Hes7 gene regulatory 

CNj ' networks. 

m ■ 

^ ; 1 Introduction 

CSI ' Periodic bodily functions are often related with oscillatory gene expression in living cells. 

In order to understand the underlying mechanism of the complex os cillatory dynamics , 
many theoretical works have been devoted over the past decades (see Klipp et al. ( 2009|) 



X 



for example). Most of the theoretical studies rely on the numerical simulations of the 
^ . detailed mathematical models of specific biological systems. However, such approaches 

i hardly provide a unified insight that is applicable to a broad class of biological networks. To 

overcome these challenges, we here consider a simplified model that is suitable for capturing 
the essence of a large class of gene regulatory networks, and develop analytic tools that are 
useful for systematically studying the existence of oscillations. 

One of t he pio neering theoretical analyses of oscillatory gene expression was presented in 



Goodwin! ( 1965[) . where the dynamical model of cyclically interconnected gene's products 



was introduced. Later, the cyclic feedback structure was found in inetabol ic and cellular 
signaling pathways as well ( Kholodenkol bOOOl : IStephanopoulos et al. . 1998 ). and it is re- 



cently considered that cyclic st ructure plays a k ey role to produce the various dynamical 
behaviors of protein levels (see Hori et al] ( 201l[ ) and referenc es therein). In fact, the ar - 



tificially constructed biological oscillator, named Repressilator (jElowitz and Leibleil [2OO0I) 



was performed with a simple cyclic interaction of repressors in Escherichia coli. These works 
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suggest that the study of cychc gene regulatory networks is the first important step toward 
the comprehensive understanding of the large-scale gene regulatory networks in nature. 

The dynamical properties of cyclic pene regulatory net works have been act ively s tudied in 



19911) to oscil- 



recent years, ranging from sta bilitv (lArcak and Son tag'. '2006', '2008'; 'ThroiJ, 
lations ()E1-Samad et al.i . .20051 : iHori et all 2011 ). In jill-Samad ct al (200j), the parameter 
conditions for the existence of oscillations were considered for Repress ilator, and t hese c on- 
ditions were generalized to a general cyclic gene regulatory networks in iHori et al.l(|201l[ ). A 
remarkable feature of these results is that the conditions are analytically obtained in terms 
of biochemical parameters despite the nonlinearity of the system. In particular, the de- 
pendence of the equilibrium point on the system's parameters is explicitly analyzed. Thus, 
one can easily see the relation between the parameters and the existence of oscillations for 
a broad class of genetic networks. As a result, no vel essential param eters that primarily 
determine the existence of oscillations were found in Hori et al] ( 201l[ ). 

In these previous works, however, the inherent time delays in transcription, translation 
and translocation process were not considered in the model. Such time delays are essential 
especially for eukaryotic cells, because mRNA and protein productions occur at different 
locations in a cell, and t he transportation of these substances results in sizable time delays 
(jChen and Aiharal 120021 ). T he existence of osc i llation s wa s studied for t he gene regulatory 
networks with time delay in Chen and Aihara (2002) and Enciso ( 2000l ). In these papers, 
however, the relation between the parameters and the existence of oscillations was not 
obtained in an explicit way. 

The objective of this paper is (i) to provide an analytic framework to study the existence 
of oscillations in cyclic gene regulatory networks with time delay, and (ii) to characterize 
the condition for the existence of oscillations with essential parameters. To this end, we 
first show that local instability of the unique equilibrium state results in oscillations of pro- 
tein concentrations based on the Poincare-Bendixson theorem for cyclic time delay systems 
(|Mallet-Paret and Selll . I 19961 ). This reduces the analysis to the local stability analysis of a 
unique equilibrium of the networked time delay systems. The main theoretical contribution 
of this paper is the derivation of graphical and analytic conditions for the existence of os- 
cillations. In particular, the dependence of the equilibrium on the system's parameters is 
explicitly considered, thus the analytic condition is obtained only in terms of biochemical 
parameters. These results ar e demonstrated with two existing gene regulatory ne t works , 



Repressilator (jElowitz and Le iblor. 2000) and the somitogenesis clock (jHirata et al.l . 120041 ). 



We also present the effect of time delays on oscillations based on the analytic conditions. 

This paper is organized as follows. In Section 2, the model of cyclic gene regulatory 
networks with time delay is introduced, and its dynamical properties are considered from 
the viewpoint of nonlinear analysis. Then, the linearized model is systematically constructed 
in Section 3. In Section 4, we derive the graphical condition for the existence of oscillations. 
Based on this result. Section 5 provides the analytic conditions. In Section 6, a toy numerical 
example is demonstrated to elucidate the developed conditions, and biological insights are 
also presented. Section 7 is devoted to the analysis of two existing biological networks. 
Finally, concluding remarks are given in Section 8. 

A part of the results in this paper was previously presented in the authors' conference pa- 
per ((Takada et all . I2OIOI ) with the omission of some details. In this version of the paper, we 
include the complete proofs of the theorems and the detailed description of the numerical 
simulations. Furthermore, the r elatio n of our developed theorem and a conflicting theo- 
retical result (Chen and Aiharal 120021 ) is clarified. Specifically, we disprove Theorem 2 in 
Chen and Aiharal (|2002h ~ and illustrate with a counter example. It is also the first time to 



present the analysis result of the somitogenesis clock in Section 7. 
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Figure 1: Cyclic gene regulatory network with time delay 



2 Modeling and Nonlinear Analysis for Gene Regula- 
tory Networks with Time Delay 

2.1 Model of cyclic gene regulatory network with time delay 

The well-known central dogma of molecular biology states that gene expression consists of 
the transcription and translation steps. During the transcription step, genes are decoded 
into molecules called messenger RNA (mRNA). Then, the information coded in mRNA is 
translated into proteins during the translation step. The rate of mRNA production is affected 
by the proteins called transcription factors, which are also created by the transcription- 
translation steps. Thus, there is an elaborate feedback mechanism to regulate protein levels 
in a cell as illustrated in Fig. [TJ This networked system is called gene regulatory network. 

In this paper, we consider the gene regulatory network, where each protein activates 
or represses another transcription in a cyclic way as depicted in Fig. [TJ In particular, 
time delays arising from transportation and intermediate chemical reactions are explicitly 
considered. The dynamics of the cyclic gene regulatory network composed of N genes is 
modeled as 

\pi{t) = -biPi{t) + Cin{t - Tr,), 

for i = 1,2,- •• ,iV, where ri,pi G M+(:= {a; G M | a; > 0}) denote the concentration of 
the i-th mRNA and its corresponding protein synthesized in the z-th gene, respectively 
( Chen and Aiharal 2002f l. For the sake of notational simplicity, we regard the subscript 



as TV and the subscript A'^ -I- 1 as 1 throughout this paper. Positive constants ai,bi,Ci and 
f3i have the following biological meanings: and bi denote the degradation rates of the i-th 
mRNA and protein, respectively: Ci and /3i denote the synthesis rates of the i-th mRNA 
and protein, respectively. A monotonic function fi : — >■ R-|_ represents either repression 
or activation of the transcription. For repression, fi{-) is defined as fi{-) := f^{-) such 
that /-(O) = 1 and /"(oo) = 0. For activation, /(•) := /+(•) such that /+(0) = and 
/^(oo) = 1. In this paper, we use the following Hill functions: 

r(p):=-l^, f+ip):=^-, (2) 
1 + p'^ l+p'^ 

where i/(> 1) is a Hill coefficient, which represents the degree of cooperative binding and de- 
termines the degree of nonlinearity of the system. The constants r^- and Tp. (z = 1, 2, • • • , N) 
represent time delays associated with the transcription and translation processes, respec- 
tively. 
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Let the following assumption be imposed throughout this paper: 
Assumption 1. 

N 



t=\ for M-)=f (■)• 

It is known that almost all solutio ns of ([I]) are observed t o asy mptotically converge to one 



of equilibria in the case of (5 > (iMallet-Paret and Sell Il996l) . Thus, it is reasonable to 



impose Assumption 1 to study the existence of oscillations, which is of our main interest in 
this paper. 

2.2 Omega-limit set of the system 

The omega-limit set of the gene regulatory network system ([Ij can be specified by using a 
Poinc are-Bendixson type theorem for time delay systems derived by iMallet-Paret and Sell 



(|l996l ) . The following proposition allows us to see that the omega- limit set of ((T|) is actually 
restricted to equilibrium points, periodic oscillations or homoclinic and heteroclinic orbits, 
and chaos is ruled out. 



Proposition 1. ([Mallet-Paret and Sdl . Il996h Consider the following syste 



ii{t)=gi{xi{t),Xi+i{t)), (i = l,2,--- ,n-l) 
in{t) = gn{x„{t), xi{t - 1)), 

where gi{-, •) (i = 1, 2, • • • , n) are nonlinear functions satisfying 



(4) 



Zi 7^ >0 and z, = <^ (5) 

av I z* e {-1, 1} if 

Let x{t) ~ [xi{t), X2{t), ■ ■ ■ ,Xnit)] G M" be a solution of ^ on some interval [tQ,oo), and 
assume that x{t) is bounded in R" as t — > oo. Then, the omega-limit set of x{t) consists of 

(a) a single non- constant periodic orbit, 

(b) equilibrium points, or 

(c) homoclinic and heteroclinic orbits. 

The dynamical model of gene regulatory networks ([T]) can be transformed to the form Q 
satisfying ([5]) by letting n = 2N and X'l cLS follows. 

\x2i^i{t) = a2i-iPN-i+i{Tt~ r]2i-i), 
1 X2iit) = a2irN-i+i{Tt - ?72i), 



for i = 1, 2, • • • , A^, where 

2N 



for i = 1 



T > T, and rii . , „, 

' IE;=2^j forz = 2,3,--- ,2iV 

with T2i-i ■■= and T2i '.— TVjv-i+i- The constants (7^ (i — 1, 2, • • • ,2N) take either 

+1 or —1, and they are defined by 

[ 1 for i = 1 

'^'■~ln-=2Pj forz = 2,3,--- ,27V, 
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where 



P2i 
P2i 



-1 := sgn 
1 



dp 



for i = 1,2, - ■ ■ ,N. 



The constant z* is then determined as z* = Y[i=iPi- The detailed proof is provided in 
Appendix |^ Note that the above transformation affects only the sign of the omega-limit 
set, thus, the omega-limit set of ri{t) andpi(i) can be specified, once that otxi{t) is obtained. 

Boundedness of x(t) is easily verified in the similar way to lEl-Samad et aL 1 (l2005l) . where 
non-delay cyclic gene regulatory netwo rks were consider ed. Existence of an equilibrium 
point and its uniqueness were proved in iHori et al.l (j201ir) in the case of no time delay, and 
time delay does not affect these properties of the equilibrium point. Hence, the following 
lemma readily follows from Proposition 1. 

Lemma 1. Consider the cyclic gene regulatory networks modeled by (Qp. Then, the protein 
levels pi{t) {i — 1,2,- •• , N) exhibit (a) non- constant periodic oscillations, (b) convergence 
to the equilibrium point, or (c) homoclinic orbits. 

Note that chaos is ruled out for the system ([1]). This lemma immediately leads to the 
following proposition, which becomes a key to deriving existence conditions of oscillations 
in Section m and [5l 

Proposition 2. Consider the cyclic gene regulatory networks modeled by (Qp. The system 
has oscillations of protein levels pi{t) {i — 1,2,- •• ,N), if the unique equilibrium point is 
locally unstable. 

In this paper, the term 'oscillations' refers to both non-constant periodic and homoclinic 
orbits. It should be noted that oscillations are periodic except for the case of homoclinic 
orbits. 



3 Linearized Model of the Cyclic Gene Regulatory Net- 
works 

We see from Lemma 1 that the local instability at a unique equilibrium point implies the 
existence of oscillations. Thus, we here consider the linearized model of gene regulatory 
networks, and derive a unified formulation that is suitable for studying large-scale gene 
regulatory networks. 

Consider a linearized system of ([T]) at the unique equilibrium point [rl,pl, ■ ■ ■ , r'^,P*N]'^ ■ 
Let fi{s) and Pi{s) denote Laplace transform of ri{t) and Pi{t), respectively. Then, Laplace 
transform of the linearized system with zero initial condition is obtained for i — 1,2, ■ ■ ■ , N 
as 



sfj(s) 






■ 




n{s) 






spi{s) 




c,e-'^^^ 








+ 






u^{s), 



where 



Ms) ■= QPi-i{s) and (i ■= fiipUi)- (7) 
Thus, the transfer function of the z-th gene from iii to pi denoted by hi{s) is computed as 



hi{s) 



{Tr,S + l)iTp^s + 1) 



(8) 



5 



where 



Ri 



Tr. 



1 



(9) 



The constant Ri represents the ratio of the geometric means of synthesis rates and degrada- 
tion rates of the i-th gene, and it is known as one of the e ssential biologica. 1 quantities which 
characterize the oscillations in gene regulatory networks (jHori et alJ . 120111 ) . The cyclic gene 
regulatory network system can be considered as the cyclic interconnection of the dynamical 
system hi{s) (i = 1, 2, • • • , N). 

In order to capture the essential dynamical properties in an analytic way, we hereafter 
simplify the model so that the kinetic parameters fli, &i, Q and Pi are homogeneous between 
genes. 

Assumption 2. We assume ai = a2 = • ■ ■ = aN{=- a), and bi = b2 = ■ ■ ■ — &Ar(=: b), i.e., 
niRNAs and proteins have common degradation rates between genes. 

With the Assumption 2, the overall system can be written by a transfer function H{s) 
defined by 



n{s):={cf>{s)e'^I-Kr\ 



1 



where 



[TrS + l){TpS + 1)' 



T 

1 -' n 



1 

6' 



(10) 



(11) 



1 ^ 



(12) 



K := 

























CnR_ 



CiRl 








(13) 



The time delays r^. and Tp. (i = 1, 2, • • • , N) of hi{s) can be different between genes, but the 
cyclic structure and the distributive property of linear systems allows us to equally distribute 
the time delays over all genes. Thus, the system 'H(s) shares the time delay t among all 
genes, where r represents the average of time delays in the network. Note that the structure 
of the matrix K is determined from the graph topology of the gene regulatory network, and 
(j){s)e'^'^ is determined from the dynamics of each gene's expression. 

In the next section, we will show that the instability of 'H(s) can be systematically checked 
using a simple graphical condition. Combined with Proposition 2, this graphical instability 
condition gives the existence condition of oscillations. 

Remark 1. Although we impose Assumption 2 to develop qualitative rather than quan- 
titative analysis framework, we can relax this assumption to some extent. In Appendix 
IB| we discuss local instability of the unique equilibrium under the relaxed assumption. In 
particular, the homogeneous case, which is considered in the following sections, is shown to 
be a worst case for local instability under parameter perturbations. Thus, the analysis of 
the simplified model can also be interpreted as the worst-case analysis. 
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Figure 2: The instability region and the eigenvalue locations of K. 



4 Graphical Condition for the Existence of Oscillations 

In this section, we present a graphical condition for the existence of oscillations. We first 
derive a necessary and sufficient graphical stability condition of H(s). Then, the existence 
condition is straightforwardly obtained, since the existence is guaranteed if H(s) is unstable 
as seen in Proposition 2. 

4.1 Local stability condition 

We first introduce an instability region of the large-scale linear system ^(s), which is char- 
acterized by the gene's dynamics /i(s)e~*'^. Let a set of complex values 17+ be defined 
as 

n+ := {A e C| 3s e C+, (f>{s)e'^ = A}, (14) 

where C+ := {s e C | Re[s] > 0}. The set i7+ is the image of the open right-half complex 
plane under the mapping (t){s)e^'^ . The instability of 'H(s) can be characterized by 57+ and 
the matrix K as follows. 

Lemma 2. Consider the system 'H{s) defined by hlO^) . Then, at least one pole of'H{s) lies 
in the open right half plane of the complex plane, if and only if 

spec{K) nn+^d). (15) 



The id ea behind this lemma is essentially the same as the one in Proposition 5.1 ofH araetal 



(l2009l ). The complete proof is presented in Appendix [C] An example of the instability region 
r2+ and the eigenvalues of K is illustrated in Fig. [2] 

The stability counterpart of Lemma 2 is characterized by the set il'^ := {A € C | Vs G 
C+, 0(s)e*'^ ^ A} with C4. := {s G C | Re[s] > 0}, which is an open complementary set of 
51+ . Then, it follows that all the eigenvalues of K lie in the stability region ff^ if and only 
if H(s) is asymptotically stable. 

Rema rk 2. Regarding the necessary and sufficient stability con dition of "H(s), Chen and Aiharal 



( 2002 ) presented a similar graphical test (see Theorem 2 in 1 Ch en and Aihara (|200^ ). The 
authors of this paper, however, have found that their graphical test is incorrect. We here 
briefly demonstrate a counterexample. More details are presented in Appendix [Dl 
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Figure 3: A counterexample to Theorem 2 in lChen and Aiharal ()2002l) 



Let iV = 3,a = 6 = 1.0000, c, = A = 1.7498, = 2.0000 and r^, = Tp, = 0.5000 {i = 
1,2,3). It follows that Ms y^ = (s + l)^e^ i?? ^ = i?| = 1.7498 and Ci = C2 = 
^3 = 0.6858. Theorem 2 in IChen and Aiharal ( 20021 ) states that the system 'H(s) is stable 



if and only if all the eigenvalues of the matrix K lie inside the region specified by an 
Archimedean spiral illustrated in Fig. 3(a) The eigenvalues of the matrix K are obtained 



3(a) 



as {1.2e^^^/^, — 1.2}, and 'H(s) would be concluded as stable from Fig. 

However, the trajectory starting near the unique equilibrium exhibits oscillations as shown 
in Fig. |3(b)[ where the initial values are set as [ri,pi, r2,p2, ''3,^3] — [0.699,1.224,0.698, 
1.226,0.697, 1.225]. In fact, we also show in Appendix ID] that the Nyquist contour of H{s) 
encloses — 1 + 7O, which impli es that H(s) is unstable (see Fig. |8]). This contradicts Theorem 



2 in IChen and Aiharal (|2002l) . It should be noted that the stability condition presented in 



this paper concludes that ■H(s) is unstable as illustrated in Fig. 3(a) 



4.2 A sufficient graphical condition 

Recall that local instability of 'H(s) is a sufficient condition for the existence of oscillations 
as shown in Proposition 2, and the instability of H{s). The following graphical condition is 
the direct consequence of Proposition 2 and Lemma 2. 

Proposition 3. Consider the cyclic gene regulatory networks modeled by (QJ). Suppose 
Assumptions 1 and 2 hold. Then, the system has oscillations of protein concentrations 
p,{t) (j = l,2,... ,iV), if 



(16) 



Proposition 3 provides a graphical condition for the existence of oscillations. Given the 
dynamics of each gene h{s)e~^'^ and the interaction matrix K, the existence of oscillations 
can be characterized by the curve defined by the inverse of h{s)e~^'^ and the eigenvalues of 
K. 

A remarkable feature of the graphical condition is that the eigenvalues are equally dis- 
tributed on a circle with the center at the origin as illustrated in Fig. |21 Specifically, the 
eigenvalues of K are written as 



Al, := Le' 



(fc = l,2. 



<N) 



(17) 
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with 



N 



(18) 



Note that L is the radius of the circle. We designate L as the average gain, since it is a 
geometric mean of the feedback gains of K in 

Thus, Proposition 3 imphes that the existence of oscihations is guaranteed if Afe e 51+ for 
some k = 1,2,- ■■ ,N. Moreover, the characteristic eigenvalue distribution allows us to 
simplify the graphical condition as follows. 

Theorem 1. Consider the cyclic gene regulatory networks modeled by (Qp. Suppose 
Assumptions 1 and 2 hold. Then, the system has oscillations of protein concentrations 
Pi{t) (i — 1,2, ■ ■ ■ , N), if \k G 51+ for some k — 1,2, ■ ■ ■ , N . Moreover, the following four 
conditions are equivalent. 

(i) Afc G 17+ for some k = 1,2, ■ ■ ■ , N . 

(ii) Ai e 51+. 

(Hi) 3wu such that |0(jwjj)e-'''^«^| < L and aig{(j){jio^)e''^*'') = t:/N. 
(iv) 3a;, such that arg ((/)(jcj,)e-'"*'^) > ■k/N and |0(jCLi,)e-'"*'^ | = L. 
where arg(-) is the argument of a complex number. 

Proof. The condition (i) is equivalent to (|16p . thus the system ^ has oscillations if (i) is 
satisfied. We hereafter show the equivalence of the four conditions. 

(i) (ii): The proof is mainly based on the fact that both |0(ja;)e^"*| and arg(0(j'a;)e^"'^) 
monotonically increase for positive w. The monotonicity is obvious from the definition (jlOp . 
Then, it is easily verified that Ai, which is the eigenvalue closest to the positive real axis, 
always goes inside the region 51+ first, since the eigenvalues of the matrix K are located on 
a circle center at the origin and radius L (see Fig. 2). 

(ii) <^ (iii): Let denote a frequency such that arg((/)(jwu)e^"'''^) — n/N. The conclusion 
immediately follows from the fact that |Ai| = L and arg(Ai) = tt/N. 

(iii) -i^ (iv): We only show (iii) (iv), since the converse can be shown in a similar manner. 
Suppose (iii) is satisfied. Let be defined such that |^(ja;*)e''"*'^ | = L. It follows that 
wjj < w*, because \(j){jojfi)e^'^>'^\ < \(j){joj^)e^'^''^\ = L and there is the gain monotonicity 
for (j){j(jj)e^^'^ as shown above. Then, the phase monotonicity implies axg{4>{juj^)e^"^^) = 



The condition (ii) in Theorem 1 greatly simplifies the graphical condition, because the 
existence of oscillations can be determined by the position of one specific eigenvalue Ai and 
the region 51+. The condition (iii) is an analytic version of the consequence (ii), though it 
is generally difficult to obtain uj^ in terms of the system's parameters. In the next section, 
the condition (iv) plays a key role to derive the analytic conditions for the existence of 
oscillations. 

5 Analytic Condition for the Existence of Oscillations 

5.1 Analytic conditions in terms of average gain 

In this section, we provide analytic existence conditions of oscillations based on the geometric 
consideration of the graphical condition. 



^/N < arg((/.(jw,)eJ''^*^). 



□ 
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We first introduce norniafized parameters of gene regulatory networks to avoid notational 
complexity, and to capture the essence of mathematical conditions. Let Ta and Tq denote 
the arithmetic and geometric means of the mRNA and protein degradation time constants, 
i.e., 

Tr+T, 



Ta:=^^^, TG:=y^. (19) 

The constants Ta and Tq have the physical dimension of time. Define the following dimen- 
sionless constants Q,ijJ and f, 

Q:^^, Co:^ojTa, f:=^. (20) 

Then, the boundary of the region il+ defined in (|14p can be written as 

0(ja;)eJ'"^ = {-Q^Cj^ + 1 + 2jw)eJ"^. (21) 

We see that the eigenvalues of K and the region are characterized in analytic form as 
([T7|) and (HH), respectively. This leads to analytic conditions for the existence of oscillations. 
We first show existence conditions in terms of the average gain L in (|17p . 



Theorem 2. Consider the gene regulatory networks modeled by (Qp. Suppose Assumptions 
1 and 2 hold. Define the two functions W{N, Q) and D{Q, L) as 



W{N, Q) =, (22) 

,s2 - +g2sin2^ 



D{Q,L) ■.= 4{l-Q^)+Q''L^ (23) 

Then, the system has oscillations of protein concentrations pi{t) (i = 1,2, • • • ,N), if one of 
the following two conditions holds □ . 

(i) L > W{N,Q), 

(a) 1<L< W{N,Q) and 



arg ( 2 - ^D{Q, L) + - 2 + ^D{Q, L) 

^ 702_2+/D(^ 
> iV - ^2^) 



Proof. It follows from Proposition 2 that there exists oscillations if the unique equilibrium 
point of ([1]) is unstable. Hence, we consider the instability condition of H(s), for which the 
simple graphical test in Lemma 2 is available. 

We first consider the case of L < 1. It should be noted that the average gain L is the 
radius of the circle where eigenvalues are located. It follows that L < 1 < |(/)(jw)e^"'^| for 
all uj. Since \(\){jLS)e^'^'^\ — 1 only when a; = 0, and a,rg{4>{juj)e^'^'^) = for a; = 0, there is 
no satisfying the condition (iv) in Theorem 1. Thus, Theorem 1 implies Ai ^ fi+, and 
T-L{s) is not unstable. 

^ This condition is necessary and sufficient for local instability of W(s), which is readily seen from the 
proof. 
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In t he case of L > W{N, Q), we readily see Ai G S1+ according to Theorem 2 in lHori et al 
(|201l[ ). In the case of 1 < L < W{N,Q), consider the condition (iv) in Theorem 1. Then, 
|<?!)(jw,)e^''^*^| = L yields 

Q^Cot + 2{2 - Q^)^l + 1 - ^0, (25) 
where :— uJt:TA- Then, is obtained as 



r-2+^DiQ,L) 
^* = Q2 . (26) 

and p4)) is derived by substituting into arg(0(jajH.)e-''^*'^) > tt/A^. □ 

The existence of oscillations can be determined by substituting the given parameters into 
Theorem 2. In particular, biological insights can be obtained by observing the relation 
between the quantities, which will be introduced in Section [6.21 

Since (I24|) has a certain monotone property in terms of L, we can simplify Theorem 2, 
and obtain the equivalent condition as follows. 

Corollary 1. Consider the cyclic gene regulatory networks modeled by (Qp. Suppose As- 
sumptions 1 and 2 hold. Define W{N, Q) and D{Q, L) as i22\) and i23\) . respectively. Then, 
the system has oscillations of protein concentrations Pi{t) (i = 1,2, ■ ■ ■ , N), if L > L, where 
L is the solution of the equation 



arg 2 - JD{Q, L) + j2x - 2 + JD{Q, L) 



vr 
iV 



Q2-2 + V^(g,L)^ 



(27) 



In particular, the solution L is uniquely determined and satisfies L € {1,W{N,Q)]. 
The proof of Corollary 1 is provided in Appendix lEl 

Re mark 3. In the case of r = 0, Theorem 2 and Corollary 1 coincide with Theorem 
2 in iHori et al. ( 2011 ). which provides an existence condition of periodic oscillations for 



non-delay case. 

5.2 Analytic conditions involving equilibrium point analysis 

In Theorem 2 and Corollary 1, the value of L depends on Q, which is determined from the 
equilibrium point p*. Since the equilibrium point depends on the parameters a,b,Ci and 
Theorem 2 and Corollary 1 require a numerical computation of the equilibrium point 
to determine Q. It is, however, desirable to explicitly take its dependence into account 
in the analytic conditions in order to gain biological insights on the relation between the 
parameters and the existence of oscillations. In this section, we restrict our attention to a 
class of the cyclic gene regulatory networks, and present analytic conditions for the existence 
of oscillations that explicitly consider the dependence of the equilibrium on the parameters. 

Specifically, we consider the case where all interactions between genes are repressive, i.e., 
fi{.) = /-(.) for alH = 1, 2, • • • , A^, and i?i = i?2 = • • • = Rn(=: R). It should be no ted 



that this class of regulatory networks includes Repressilator ( Elowitz and Leibler . 20001 ). 



We first see that the equilibrium p* does not change by ti me delay, and p* = p^ = ■ ■ ■ = 
p*^{^ p*) and d ^ Q* ^ ... = Q{= (*) hold as shown in iHori et ail (|2010[ ). Using this 
property, we have the following relation between L and R. 
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Lemma 3. Consider the gene regulatory networks modeled by (QJ). Suppose fi{-) = 
/"(•) {i = 1,2,--- ,N), Ri = i?2 = ••• = RNi—- R), and Assumptions 1 and 2 hold. 
Then L < v, and the following equality holds. 

The proof of this lemma is provided in Appendix IF] 

Lemma 3 shows a direct relation between the average gain L and the biological parameters, 
R and v. Then, this lemma leads to the following existence condition of oscillations, which 
is explicitly written in terms of the biological parameters. 

Theorem 3. Consider the gene regulatory networks modeled by (Qp. Suppose fi{p) = 
f~{p) {i = 1,2,--- ,N), Ri — R2 — ■■■ — i?Ar(= R), and Assumptions 1 and 2 hold. 
Define W{N,Q), D{Q, L) and L as if^j . ifI3) and (F^ , respectively. Then, the system has 
oscillations of protein concentrations Pi(t) (i = 1, 2, • • • , N) if both u > L and R > R hold, 
where 

«^-(^)"^^ (-) 

Proof. We derive an equivalent condition to Corollary 1. Observe that R^ in p8| is 
monotonically increasing for L{< v). Thus, if the condition L < L in Corollary 1 is satisfied, 
R < R follows, where R^ is defined by ([29l) . We see from Lemma 3 that > L is also satisfied, 
because v > L. On the other hand, if > Z and R > R are satisfied, we have L > L because 
of the monotonicity of ([25]) . 

Since the conditions v > L and R > R are equivalent to those of Corollary 1, we can 
conclude that there exists oscillations if these conditions are satisfied. □ 

Theorem 2 provides a condition for the existence of oscillations in terms of biological 
parameters z/, i? and R{i>, L{N,Q,f)) without any information about the equilibrium point 
p*. This is contrast with the conditions in Theorem 2 and Corollary 1, where the dependence 
of L on the equilibrium point p* is not explicitly obtained. Therefore, we can conclude from 
Theorem 2 that the following five dimensionless parameters characterize the existence of 
oscillations: the number of genes (N), time delay normalized by the arithmetic mean of the 
lifetime (f), the Hill coefficient (v), the ratio between the geometric mean of degradation 
rate and production rates (i?), and the ratio between the geometric and arithmetic means 
of degradation rates (Q). 

Remark 4. We can see that R and L are monotone with respect to the system's parame- 
ters. This observation leads to the conclusion that the system tends to have oscillations by 
increasing any of the five essential quantities, N, f , v, R and Q. We will study more details 
in Section 15^ 

6 Numerical Example and Biological Insight 

6.1 A cyclic gene regulatory network with N = 7 

We here confirm the theoretical results provided in Theorems 1, 2 and 3 with illustrative 
numerical examples. 

Consider the cyclic gene regulatory networks composed of A'^ = 7 genes. Assume that 
a = L2,6 = 4.8, ci = C3 = ce = C7 = 1.92, ca = C4 = C5 = 3.84, ^1 = /Jg = = /?? = 
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(a) Case: r = 1. The protein concentra- (b) Case: r = 0. The protein concentrations 
tions exhibit oscillations. converge to the equilibrium. 



Figure 4: Time plot of protein concentrations. 



4.32,^2 = /34 = /35 = 2.16, and let the Hill function be defined as f^{p) = 1/(1 with 
v — 2.6 for i = 1, 2, • • • , A^. Then, Q and R{:— Ri — R2 = ■ ■ ■ = R7) are obtained as 
Q ^ 0.800 and R = 1.200 from the definition ^ and (H]), respectively. The value of L 
in (|17p can be computed as L = 1.048. Note that Q in L involves computation of the 
uni que equilibrium of the system, but it can be efficiently done with the bisection algorithm 
(see Hori et al. ( 2010[ ) for details). In the following, we will see the effect of time delay by 
comparing a genetic regulatory network with and without time delay. 

We first apply the graphical existence condition in Theorem 1 . Figure [2] illustrates the 
instability region $7+ and the eigenvalue distribution oi K ior f ~ and f ~ 1.00, where 
the time delays of the reactions are set as t^j = r^g = r^^ — — 0.31, — — Trg = 
0.26, — Tp^ — Tp^ — Tp^ — 0.21, = — Tpg = 0.26, thus the average of the time 
delay is r = 0.52. In the case of f = 1.00, the boundary of the instability region r2-|_ 
is given by (/)(jw)e^"'^ in Fig. [2] Thus, Theorem 1 implies the existence of oscillations, 
because two eigenvalues of K belong to the region Q^. In the case of f = 0, the boundary 
of the instability region il^ is 4'{j(jj) in Fig. O We can see that all eigenvalues of K are 
located outside the region f7+ when f = 0. Thus, it is concluded from Theorem 1 that 
a unique equilibrium point of the system is locally asymptotically stable, and the protein 
concentrations do not exhibit oscillations when they are perturbed around the equilibrium 
point. Note that this result does not imply non-existence of oscillations, since Theorem 1 is 
a sufficient condition for the existence of oscillations. 

The same conclusion follows from the analytic conditions in Theorem 2. We can see from 
Theorem 2 that there exist oscillations when f ~ 1.00, because L ~ 1.048 > L = 1.031, 
where L is computed by ([?7|) . On the other hand, L = 1.048 < L = 1.072 in the case of 
f = 0. Since the condition in Theorem 2 is equivalent to that of Theorem 1, we can conclude 
that the equilibrium point is locally asymptotically stable. 

Theorem 2 and Corollary 1 required the value of equilibrium point to compute L. In 
contrast. Theorem 3 does not require the computation of the equilibrium. For given param- 
eters, L and R can be determined from ((27|) and (|29p . respectively. Specifically, L = 1.031 
and R = 1.187 when f = 1.00. Computing R from ^ yields R = 1.200. Therefore, both 
of 1/ = 2.6 > Z = 1.031 and R = 1.200 > ^ = 1.187 in Theorem 3 are satisfied, and 
the existence of oscillations is concluded. On the other hand, R — 1.218 when f = 0. 
Thus, the conditions in Theorem 3 do not hold, because R = 1.200 < R — 1.218 despite 
v = 2.6 > L — 1.072. It should be noted that this result, in turn, implies that there exists 
oscillations even for f = 0, if the parameters a,b,Ci and /3i are set so that R > 1.218 and 
L < 2.6. 



13 










1 




10 




Existence of oscillations 



3 4 
V 




— 2=0.1 

— 2=0-5 

— e=i.o 



Existence of oscillations 



1_L 



3 4 
V 



(a) f = 0, 1 and 10 with = 3 and Q = 1 (b) Q = 0.1, 0.5 and 1.0 with = 3 and 

f = 1 




(c) TV = 3, 5 and 11 with 



Figure 5: Parameter regions {v, R) for the existence of oscillations. 

In fact, numerical simulations shown in Fig. |4]show oscillations and convergence to the 
equilibrium of protein concentrations for the time delay and non-delay case, respectively. 
We have seen that the existence of oscillations is more probable when the time delay is 
large. In what follows, we will see that this is indeed the case in general. 



6.2 Biological insight: relation between parameters and existence 
of oscillations 

As we have seen in Theorem 3, the existence of oscillations in cyclic gene regulatory networks 
with time delay can be characterized by the five dimensionless parameters N, Q, f, v and R. 
The parameter regions that guarantee the existence of oscillations can be drawn as shown 
in Fig. [5] based on the analytic conditions given in Theorem 3. From these figures, we can 
readily conclude that the system tends to have oscillations as ly and R get larger. In addition, 
the larger the parameters f, Q and N are, the more probable the existence of oscillations 
becomes because of each figure in Fig. [5] 

An advantage of Theorem 3 is that we can confirm that these observations are true in 
general because the conditions are written in an analytic form in terms of the given biological 
parameters. Therefore, we conclude that the gene regulatory networks, in general, tends to 
have oscillations by letting the five essential parameters TV, Q, f , v and R larger (see also 
Remark 4) H. 



^The constant R in Theorem 3 is not monotone decreasing for all u (> 1), but it becomes a decreasing 
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Figure 6: Parameter regions (a,7~^) for the existence of osciUations in Repressilator. Both 
axes are in common log scale. Oscillations are more probable as the time delay becomes 
large. 



7 Applications 

In this section, we will apply our results to two existing biological systems, and see how our 
results work in analyzing the effect of time delay on the existence of oscillations. 



7.1 Repressilator 



Repre ssilator is one of the pioneering synthetic gene regulatory networks created bv lElowitz and Leibler 
(|2000l ). This artificial cyclic gene regulatory network is composed of three repressor genes, 
each of which represses another gene and forms cyclic reaction structure shown in Fig. 
[TJ In lElowitz and Leibleil (|2000f) . Repressilator was implemented in Escherichia Coli, and 
oscillations of protein concentrations were observed in vitro. 

The dynamical model of Repressilator can be written as 



rAt) = -rJt) H — 



(30) 



for i = 1,2, 3, where 7 denotes the ratio of the protein degradation rate to the mRNA degra- 
datio n rate, and the constant ao represents leakiness of the promoter (lElowitz and Leibler 
2OOOI) . Note t hat time delays are not co nsidered, i.e., Tr- — Tp. — for i ~ 1,2,3, in the 



original paper (jElowitz and Leibleil . I2OOOI) . We remark that the recently proposed technique 



function for 1 < < 8, which is the region of our interest. 
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(Ugander, 2008[) could enable u s to engineer the ti me delay, and it would contribute to ob- 
taining a desired dynamics (see Orosz et al. ( 2010f) for example). Hence, the time delays in 
the above model should account for such engineered delay as well as fast dynamics omitted 
in the modeling process. It can be seen that the model (pO]) is equivalent to ^ by rescaling 
the parameters when ao — 0. We notice that Proposition 2 holds, even when ao ^ 0, and 
Theorem 2 and Corollary 1 can be applied to analyze the existence of oscillations. 

Let us first consider the case where no time delay appears in dynamics, i.e., r^. = Tn.- 



' I 1, 2, • • • , TV), which is the original model of Repressilat or presented inlElowitz and Leibler 

Following the numerical simulations conducted in lElowitz and Leibleil (|200oir we 
set a = 624, ao = 0.0866, 0.200 and ly = 2.0. Then, L and L in Corollary 1 can be 
computed as L = 1.833 and L — 1.519, respectively. Thus, we conclude t he existence of oscil- 
l ation s from Corollary 1, which is consistent with the simulation result in lElowitz and Leibler 
(l2000l ). 

Next, we investigate the effect of time delay on the existence of oscillations, and show 
that time delay increases robustness of Repressilator. Here we numerically examined the 
existence conditions in Theorem 2 for various time delays and parameters. The result is 
shown in Fig. [SI where the parameter regions for the existence of oscillations are illustrated. 
Note that only the normalized time delay f rather than each time delay itself affects the 
existence of oscillations as seen in Section 16.21 Also note that the p a ramet er region for 
f = in Fig. [5] coincides with that in Fig. 1 (b) in Elowitz and Leibler ( 2000l ). We can see 
from Fig. [5] that the regions for the existence of oscillations get larger as f become larger. 
This implies that one could make robust oscillator by inserting time delay. Moreover, the 
parameter region is not sensitive to a little change of ao and v when f is large. 



7.2 Somitogenesis clock 



Somitogenesis is the process by which the somites of living organisms are created. Biological 
experiments as well as theoretical studies showed that the timing of the somite segmentation 
is regulated by an oscillatory expression of Hes7 gene fsee lHirata et all (|2004l ): lLewisl (|2003l) 
and the references therein). In particular, it was shown by a biological experiment that 
oscillations produced by nega tive self-feedback o f Hes7 play a crucial role in controlling 
the somitegenesis oscillations ( Hirata et al. . 2004 ). In this section, we focus on the Hes7 



regulatory networ k, and see the validi ty of our theorems by comparing with the experimental 
data presented in Hirata et al.l ( 20041) . In addition, we provide some insights obtained from 
the developed theorems. 



Following lHirata et al.l (|2004[ ). we consider the following dynamical model of the regulatory 
network of the Hes7 protein. 



r(t) = -ar{t) + 



/3 



bp{t) + cr{t ~ Tr). 



(31) 



This model is equivalent to the model (HJ setting = 1 and v = 2. Here, the mRNA and 
protein degradation rates a and h are defined as 



log 2 ^ log 2 



(32) 



where and tp denote mRNA and protein half-l i fe tim e. We employ the parameter 
values for wild- type Hes7 provided in iHirata et al.l (|2004l ): tp = 20 [min],^^ = 3 [min]. 



a = 0.231 [min" 
7 [min]. 



^],b = 0.0347 [min-i], c = 4.5 [min"!],/? = 33 [min^^], Tp 
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Table 1: Values of the essential biological parameters in somitogenesis oscillator 



Parameter 


Before mutation 


After mutation 


N 


1 


1 


Q 


0.674 


0.575 


f 


2.23 


1.55 


V 


2 


2 


R 


21.5 


26.4 


R 


6.99 




L 


1.97 


1.97 


L 


1.85 


2.39 



25 




mRNA half-life, min 



Figure 7: Parameter region for the existence of oscillations in somitogenesis oscillator. The 
existence of oscillations is mo re probable when m RNA half-life time is small, which is con- 



sistent with the hypothesis in iHirata et al.l (|2004l) 



In lHirata et al.l (|2004[ ). a point mutation in the gene was introduced, and mice expressing 
mutant Hes7 were generated. The protein half-life of one of the Hes7 mutants was identified 
as almost tp — 30 minutes, which is longer than that of the wild-type Hes7, which is 
tp = 20 minutes. As a result, the protein degradation rate of the mutant Hes7 changed 
to 6 = 0.0231 [min-i]. 

Numerical simulations of the model ([5T|) revealed that the protein of the wil d-type Hes7 



show s oscillations, but that of the mutant Hes7 converges to a stable equilibrium (jHirata et al 



2004). In addition, the experimental result was consistent with the numerical simulations 



( Hirata et al.l . l20O 



We here present that Corollary 1 and Theorem 3 can explain these observations. First, we 
compute the values of the essential biological parameters, and obtain Table 1. We see from 
Table 1 that there exist oscillations before the mutation, because L — 1.97 > L — 1.85 in 
Corollary 1, and equivalently R = 21.5 > R = 6.99 in Theorem 3. On the other hand, the 
equilibrium point can be found to be locally stable after the mutation, because L = 1.97 < 
L — 2.39 in Corollary 1, and L = 2.39 > = 2 in Theorem 3. We see that these results 
agree with the existing experimental work addressed above. 

that short half-life time t„ 



It was concluded in iHirata et al.l (|2004l ) that short half-life time tp of Hes7 protein is a 
key to the oscillations, though theoretical analysis was not performed. This hypothesis can 
be theoretically verified by using the presented theorems. Using Corollary 1, we can obtain 
the parameter region for the existence of oscillations in terms of half-life time of mRNA 
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tr and protein tp in Fig. [71 We see that robust oscillations are guaranteed if the protein 
half-life time is shortened. For example, when mRNA half- life time is = 3 [min], there 
exists oscillations for 0.1 [min] < ip < 22 [min]. 

8 Conclusion 

In this paper, we have considered the conditions for the existence of oscillations in cyclic 
gene regulatory networks with time delay. Based on the unified analysis framework of 
gene regulatory network, we have first derived the graphical condition. Then, the geometric 
consideration has led to the analytic conditions. In particular, the relation between the equi- 
librium point and the parameters are explicitly considered, thus the condition is explicitly 
written in terms of biochemical parameters. Thus, biological insights can be easily obtained, 
and the relation between the parameter and the existence of oscillations has been revealed. 
Finally, we have confirmed that the developed theorems can be helpful to determined the 
existence of oscillations in two existing biological networks. 
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A Transformation of the System 

We here show that the gene regulatory network system defined by (H)) can be obtained from 
(|4|) by the transformation ([6]). 
We take a time derivative oi Xi{t) {i — 1,2, ■ ■ ■ , 2N), and substitute ([T|). 



X2i-i{t) = a2i-iTpN-i+i{Tt - ■r]2i-i) 

= 0-2i-lT{'~bN-i+lPN-i+l{Tt — 772i-l) + 



CN-i+irN-i+l{Tt ~ r]2i-l — Tr 
= —bN-i+lTX2i-l{t) + 



)} 



CN-i+iTa2i-irN-i+i{Tt ~ r]2i) 

= —bN-i+\Tx2i-l{t) + CN-i+lTp2iX2i{t), 



(33) 



for i = 1,2, • • • ,iV. 



X2i{t) = (J2iTrN-i+i{Tt - r]2i) 



= o-2iT{—aN-i+irN-i+i{Tt — r]2i) + 

l3N-i+lfN-i+l{pN-i{Tt - r]2i - T-p„_ J)} 



= -aN-i+iTx2i{t) + 

l3N~i+iTa2ifN~i+i{pN-i{Tt — r?2i+i)) 



= -aN~i+iTx2i{t) + 

PN~i + \T(J2ifN-i + \ {cy2i + \X2i + \{t)), 



(34) 
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for j = 1,2, ■ ■ • , - 1, and 



i2]v(t) = a2NTri{Tt - ?72jv) 

= a-2NT{—airi{Tt — r]2N + 

Pl<J2N fl(pN{Tt — 772JV — Tp„))} 

= -axTx2N{t) + l3iTcj2Nh{PN{Tt - T)) 
= -aiTx2N{t) + l3iTa2Nfi{xi{t - 1)). 



(35) 



We see that (gSl), (EH) and (|3S]) are of the form (1). Also we can verify that (021), (EH) and 
([551) satisfy (5) as follows. It holds that 



gg2i-i(x2i-i,a:2i) 



CN-i+lTp2i, 



dx2i+l 
dg2N{x2N,Xi) 



PN-i+lTa2iCf2i+l 



dfN- 



i+1 



dp 



dxi 



= j3iTa2N 



dp 



(36) 
(37) 
(38) 



It is clear from the definition that (|36| is positive. We can see that (|37)) and ((38)) are also 
positive, because it follows that 



o-2icr2i+iSgn 



dfN- 



dp 



> 0, cr2AfSgn 



dh 



dp 



-z . 



Note that the sign of xi is the same as that of CTi, and /i(-) is a monotonic function de- 
fined on positive orthant. The relation ([5]) can be alternatively confirmed by the following 
calculation. 



dg2i-^\(x2i^\,X2i) _ drpf-i+i dg2i-i 



dX2i 

dg2i{x2i,X2i+l 
dX2i+l 



dX2i dVN-i+l 
CN~i+lTp2i, 

dpN-i dg2i 



(J2iCN-i+lT(J2i-l 



dX2i+l dpN^i 
~ 0-2i+ll3N~i+lTa2'. 



dp 



dg2N{x2N,Xl) dpN dg2N r, rr. dfi 

— 7^ = O-lPil (J2N-r- 

OX\ opN dp 



dxi 



PiTa2N^r'- 
dp 



where the right-hand sides coincide with those of ([55]). (1571) and ((551) . respectively. 



B Discussions on Robustness of Local Instability 

In this section, we relax Assumption 2 and discuss local instability of the unique equilibrium 
of (|T|) (see also Remark 1). Specifically, we consider robust instability of the linearized 
system under parameter perturbations, which are often the case in biochemical systems. 
Let V denote a set of parameters defined by 

V := {{a^,bi,Ci,l3^,Tr.,Tp^,Ci) = 1, 2, • • • , A^) | Oj < < al,bi <bi < b^.c^ < Ci < cl, 

ft <A< A,^<T., <—^,Tp^<Tp^ <T^,Q< \Q\ < Ci (* - 1,2,--- ,N)}, (39) 

where the symbols with a underline and an overline, i.e., ai and al etc., represent given upper 
and lower bounds of the parameters, respectively. The perturbation of the linearized gain 
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Ci (i = 1, 2, • • • , N) accounts for both the uncertainty of equihbrium due to the perturbations 
of {tti, bi,Ci, Pi) and that of the Hih coefficient v. We define the Unearized system by H(s) := 
(/ - H{s)M)-'^, where H{s) := diag(/ii(s), /i2(s), • • • ,hN{s)) and 



M 






• 


• 


Ci " 


C2 


• 


• 








C3 • 


• 








• 








(40) 



In iHori et alj (|201lh . it was shown for t^- — Tp. — Q that the worst-case perturbation 
for local instability is given by (a^, bi, Ci, Ci) = (a7, ^i, Pi,SiQ), which is the upper and 
lower extremum of the parameter set. We can easily obtain a similar result for r^. ^ and 

Tp. ^ 0. 

Proposition 4. Consider the cyclic gene regulatory networks modeled by {Ip. Suppose 
Assumption 1 holds. Then, the linearized system 'H{s) is unstable for all parameters in V , 
if and only if 'H{s) with 



fli = flj, bi = bi, a = Ci,f3i^ /3i, 



^Pi 1 Ci — ^iCi 



(41) 



is unstable. 



Proo f. The idea of the proof is essentially the same as the one presented in iHori et al 
pOllh . but we here provide the detailed proof for completeness. 

The poles of 'H(s) are given by the roots of 



N 



\I-His)M\^l + l[h,{s)\C, 



i=l 



where 7(5) := {s + bi)Ylf^^{l/hi{s)). We first fix the parameters ai{i — 2,3, •• • ,N) and 
bi, Ci,l3i, Tr-, Tp. {i = 1,2, ■ ■ ■ , N). Let ap{— ai) and Wp(= uj) satisfy P^ . which implies 

v*p = -^Jal + uj'^p\j{jujp)\ and Z{ap + jujp) + ujp{Tr, + Tp^) = n - Zj{jujp), (43) 



ci/3i 



TV 



(s + ai)7(s) + []|C,| =0, (42) 



i=l 



where v* := YiiLi K*! critical gain for destabilizing the closed-loop system. 

Here, we perturb Op to a^{> Op). Then, it follows that 

Z(a^ +jOJp) +UJp{Tr^ +rp„) < Z(a^ + jWp) +UJp{Tri +Tp^) =71 - Zj{jUp). 



(44) 



Moreover, there exists w^(> Wp) such that Z{ai, + jcui,) +uj^{Tr-^ +'^pjv) = ^ ^lU^v): since 
7(jcj) is an increasing function in terms of lo. Then, the critical gain w* that destabilizes 
the closed-loop system is given by 



1 



■Vai+^\l{j^u)\ > Wp, 



(45) 



because > Wp and |7(ja;)| is an increasing function for cj > 0. 

Thus, the critical gain, which destabilizes the closed-loop system H{s), monotonically 
increases as oi increases. This implies that H(s) is unstable for all ai € [ai,^]"], if and only 
if the system with oi = is unstable. The same conclusion follows for (i = 2, 3, • • • , N) 
and &i(« = 1, 2, • • • , A^) by the same proof. 
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Regarding a and /3i {i — 1,2, N), we see that only the gain condition in (^5]). which 
is the first equahty, depends on Ci and f5i. Thus, we can immediately conclude that Ci — Ci 
and /?i = f3i are the worst-case parameters for instability of 'H(s). 

We can apply the same discussion to the delays. Suppose Tp{— r^J satisfy We define 
a perturbed delay as t^{> Tp). Then, Z(ai + jtOp) + ujpir^ + Tp^) > tt — ^{jujp), and we 
see that there exists lu^, < ujp such that Z{ai + juji,) + Wj,(t,^ + Tp„) = vr — ^{jijj^). Thus, 
the critical gain for instability decreases as r^j increases. This, in turn, means that r^^ is 
the worst-case delay for robust instability. Applying the same proof for r^. (i = 2, 3, • • • , A'') 
and Tp. (j = 1, 2, • • • , N), we have the conclusion. □ 

This proposition means that we can guarantee the existence of periodic oscillations for all 
parameters in V only from the local instability analysis for the extreme parameters. 

Let a := max^ Hi, b := max^ bi, c := min^ Ci, jS := min; Pi, := min^ and Tp := min^ Tp. . 
Proposition 4 immediately leads to the following sufficient condition for local instability. 

Corollary 2. Consider the cyclic gene regulatory networks modeled by ([IJ). Suppose 
Assumption 1 holds. Then, the linearized system 'H(s) is unstable for all parameters in V , 
if'H{s), which is defined in U0\). is unstable for 

a ^a,b = b,Ci = c,(3i ^ P, Tr^ = tv, Tp^ = Tp, and Q = SiQ. (46) 

This corollary means that the equilibrium point is locally unstable for all V, if the homo- 
geneous system H{s) defined by (fTUl) is unstable for the extreme parameters (^5)) . Note that 
local instability of the equilibrium is a sufficient condition for the existence of oscillations as 
shown in Proposition 2. Although Assumption 2 is made in order primarily to simplify the 
model and gain qualitative insights, Corollary 2 implies that the analysis of the simplified 
model also provides conditions for the robustness of the existence oscillations. 



C Proof of Lemma 2 

The poles of 'H(s) are obtained by solving 



\(f>is)e'^I -K\^ \<i){s)e'^I - A| = 0, 



(47) 



where A := diag(Ai, A2, • • • ,\n) G C . It should be noted that 'H{s) is a retarded time 
delay system, thus the dominant pole is located at the rightmost position in the complex 
plane. Thus, it follows from the definition of f2+ that T-L{s) has at least one pole in C+ if 
and only if spec(i^) n ri+ 7^ 0. □ 



D 



A Cou nterexample to Theorem 2 in lChen and Aihara 
( I2QQ2I ) 



In this section, we first examine the stability of T-L{s) by numerically comp uting the Nyquist 
plot a nd the poles of T-L{s). We then point out the errors of the proof in IChen and Aihara 
(l2002l ). 

Let the parameters of ([T|) be set as the ones in Remark 2. It follows that (j>{s)e^'^ = 
(s + 1)^6", Rl = Rl=: Rl^ 1.7498 and Ci = C2 = C3 = 0-6858. In what fo llows, we show 



that H(s) is actually unstable, though Theorem 2 in lChen and Aiharal (|2Q02I ) concludes that 
it is stable as explained in Remark 2. 
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Figure 8: Nyquist contour of the system. The contour encloses -1+jO. 
lOOr ■ : 
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Figure 9: Roots of the characteristic equation in terms of the eigenvalue ^/l.2e^^^^. 



Figure [8] illustrates the Nyquist plot of the loop transfer function of 'H(s), which is defined 



by 



n 



Ci/3i 



-3s 



(s + l)6 



We see that the Nyquist contour encloses —1 + jO. Since the open loop system is stable, the 
Nyq uist contour in Fig . [51 im plies that the system is unstable, which contradicts Theorem 
2 in IChen and Aihara ( 2002 ). In fact, a pair of the poles of 7^(s) is found in the open 
right-half complex plane at 0.0212 ± 0.3 6347 by numerical computation (see Fig. [SJ. These 
observations imply that Theorem 2 in Chen and Aihara ( 20021 ) is not the necessary and 
sufficient stability condition for H{s). 

We hereafter clarify the errors of their m athematical proof. There are essentially two 
er rors in the m athematical proof provided in Chen and Aihara (12002 ). First, Theorem 2.6 
m lBelaiij(ll993l) . which is use d in the proof in lChen and Aiharal (|2002|). is incorrect. S econd, 



Theorem 2.6 in iBelaiij (jl993l) was applied in a wrong way in lChen and Aiharal (|2002f) 



In th e remaining of this section, we use the notations defined in lBelaiij (|1993l ): IChen and Aihara 
(|2002l) for the sake of easy comprehension and comparison. Our first claim is that Theorem 
2.6 in iBelaiii (|l993f) is not the necessary and sufficient, but a sufficient condition. Using the 
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R = ppJ^-i) 



notations in Belair (Il993l) . we see that 



<j) = 6 + ST + 27rj 

R = pe('-i)^ (a) 
(P = 9 + Re^^-''^^T sin 6 (b). 

where p > and < < 27r[l. Note that the equations (a) and (b) are the same as (2.7a) and 
(2.8) in lBelaiii (|l993t ). It was concluded that ah roots A of the above equation have negative 
real parts if and only if (i?, </>) G {{R, </)) e [0, oo) x [0, 27r) | (a) and (b) only if r < 1}. Then, 
the stability region was considered by specifying (i?, </>) that belongs to the above set. It 
follows that 



{{R, (f)) £ [0, oo) X [0, 27r) | (a) and (b) only if r < 1} 
3 {(i?, 0) G [0, oo) X [0, 27r) | (b) only if r < 1} 
= {{R, (f) e [0, oo) X [0, 27r) | (b) for some r > 1}, 



(48) 
(49) 

(50) 



where I ■} de notes a complementary set. Then, the set ([SO]) was specified in Theorem 2.6 in 
Belaiil ()l993l ) as 



{(i?,0) e [0,oo) X [0,27r) I (b) for i 



>1}, 



(i?, <?!)) G [0, oo) X [0, 27r) ] <^<^ + i?r or (I)>^-Rt 
(i?, G [0, oo) X [0, 27r) 1 <^>| + i?r and (P<^-Rr 



(51) 



In lBelaiii (|l993f) . however, was not derived as a subset but as an equivalent set of pS)) . 
Therefore, it was concluded that (ISTI) provides the stability region where all the eigenvalues 
of the system are located, if and only if the system is asymptotically stable. This conclusion 
is, however, incorrect, because (j49|) is actually a subset of (|48l) . 

Instead, we see that (|5ip is the region where all the eigenvalues of the system are l ocated, 
if the system is asymptotically stable. Therefore, we claim that Theorem 2.6 in I Belair 



(|l993f ) provides only a sufficient condition for stability. 

The other error of Theorem 2 in Chen and Aiharal (|2002[ ) stems from the mis-application of 
Theorem 2.6 inlBelaiii (|l993l ). The boundary of the stability region provided in Theorem 2.6 



m 
in 



Belaii 



Belaii (1993) is the Archimedean spiral sta rting from the ori gin. Applying Theorem 2.6 



( 1993 ) to the equation (11) in i Chen and Aiharal ( 2002 ). we see that the boundary 



of the stability region is given by 



R 



29 



fcr 

377-20 



kr 



for - < 
for TT < < 



< TT 



(52) 



where the constants k and r are defined as in IChen and Aiharal (|2002l) . The tuple {R,e) 
defines the distance and the angle of the boundary measured from the origin. Consequently, 
the arc drawn by (j52p becomes the well-known Archimedean spiral. 



^Equation (2.4) in lBelairl l ll993l) is typo. It should be corrected as A = —1 + be 
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In Theorem 2 of IChen and Aiharal (j2002l ). however, the left-hand side of ([5^ . i?, was 
shifted by one, and the equation of the boundary was given by 



R-1 



for - < 6* < TT 

fcr 2 

37r-26' ^ ^ Stt 

^ for.<^<-. 



(53) 



T hen, -R and 9 are measu red from the origin and 1 + jO, respectively (see also Fig. 2 



Chen and Aiharal (|2002h 'l. It is clear that the bou ndary obtained in this way does not 



coincide with the one in Theorem 2.6 in iBelaii ( 19931 ) 



E Proof of Corollary 1 

We observe that the left-hand side of (f24|) is the monotonically increasing function of L, 
and the right-hand side is the monotonically decreasing function of L. Moreover, we can see 
that the inequality (f24| is not satisfied when L < 1, but it is satisfied when L > W{N, Q). 
Therefore, we have a critical value L at which the left-hand side and the right-hand side 
of p4| take the same value, and the inequality p4)) is satisfied if and only if i > Z. It is 
clear from the above argument that L is given as the unique solution of (|27|). and L satisfies 
1<L< W{N,Q). 

F Proof of Lemma 3 

It follows that 



L = -RX = -R' ^ ^ , (54) 



where the second equality follows from the definition of (. According to Hori et al. ( 20101 ). 
we have 



n2 

p* = — . (55) 



Then, it follows from 1^ and that 



L^^. (56) 



This implies L < v. In addition, it follows that 



C^^^,{R'-pn (57) 



(see iHori et~al] (|2010h for the details). Muhiplicating i?^ to ([ST]) , we have 



Thus, we can eliminate p* from ([56)) by using (|58p . and obtain (28) 
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